package crack;

import inf.math.LinAlg;

import java.io.IOException;
import java.util.ArrayList;
import java.util.Collections;

import math2.ConstantFunctionRnToR;
import math2.FunctionRnToR;
import math2.MathUtilities;
import fem2.AbstractStructuralCrackDemo;
import fem2.Constraint;
import fem2.Debugger;
import fem2.DistributedLoad;
import fem2.Edge;
import fem2.Element;
import fem2.Face;
import fem2.Load;
import fem2.MaterialModel;
import fem2.Mesh;
import fem2.MeshGenerator;
import fem2.MeshUtilities;
import fem2.Model;
import fem2.Node;
import fem2.analysis.Analysis;
import fem2.enu.EchoLevelType;
import fem2.enu.SortingCriteria;
import fem2.enu.State;
import fem2.material.StVenantKirchhoffMM;
import fem2.material.cohesive_law.CohesiveLaw;
import fem2.mesh_and_topology.SpatialCriteria;
import fem2.observer.StressObserver;
import fem2.observer.TimeStepObserver;
import fem2.observer.XYGraphicalObserver;
import fem2.observer.XYSourceObserver;
import fem2.pre_and_post.GidMeshGenerator;
import fem2.strategies.Strategy;

public abstract class AbstractFourPointTestUsingSDAElementDemo extends AbstractStructuralCrackDemo {

	protected double loadFactor = 0.13;

	protected double hp = 20;
	protected double wp = 40;

	protected boolean enableFixOnePointOnLeft = false;
	protected boolean enableSlideOnePointOnRight = true;

	public AbstractFourPointTestUsingSDAElementDemo() {
		projectDir = "/home/hbui/kratos_janosch";
		projectName = "four-point-test.gid";
	}

	@Override
	public Mesh createMesh() {
		String fn = projectDir + '/' + projectName + '/' + meshFileName;
		MeshGenerator mg = new GidMeshGenerator(fn);
		Mesh mesh = mg.getMesh(2);
		return mesh;
	}

	@Override
	public Model createConditions(Model m) {
		Mesh mesh = m.getMesh();

		Constraint cxy = new Constraint(false, false);
		Constraint cy = new Constraint(true, false);

		final double d2 = 600;
		SpatialCriteria steel_cap1 = new SpatialCriteria() {

			@Override
			public boolean isSatisfied(double... x) {
				boolean cy = Math.abs(x[1] - (-hp)) < MathUtilities.tol8;
				boolean cx = (x[0] >= d2 - wp / 2) && (x[0] <= d2 + wp / 2);
				return cx && cy;
			}
		};

		if (enableFixOnePointOnLeft) {
			MeshUtilities.seekNodeNearest(mesh, steel_cap1, d2, 0).setConstraint(cxy);
		} else {
			mesh.setConstraint(cxy, MeshUtilities.seekNodes(mesh, steel_cap1));
		}

		final double d5 = 1119;
		SpatialCriteria steel_cap2 = new SpatialCriteria() {

			@Override
			public boolean isSatisfied(double... x) {
				boolean cy = Math.abs(x[1] - (-hp)) < MathUtilities.tol8;
				boolean cx = (x[0] >= d5 - wp / 2) && (x[0] <= d5 + wp / 2);
				return cx && cy;
			}
		};

		if (enableSlideOnePointOnRight) {
			MeshUtilities.seekNodeNearest(mesh, steel_cap2, d5, 0).setConstraint(cy);
		} else {
			mesh.setConstraint(cy, MeshUtilities.seekNodes(mesh, steel_cap2));
		}

		final double d1 = 203;
		final double h = 306;
		SpatialCriteria steel_cap4 = new SpatialCriteria() {

			@Override
			public boolean isSatisfied(double... x) {
				boolean cy = Math.abs(x[1] - (h + hp)) < MathUtilities.tol8;
				boolean cx = (x[0] >= d1 - wp / 2) && (x[0] <= d1 + wp / 2);
				return cx && cy;
			}
		};

		final double d4 = 722;
		SpatialCriteria steel_cap3 = new SpatialCriteria() {

			@Override
			public boolean isSatisfied(double... x) {
				boolean cy = Math.abs(x[1] - (h + hp)) < MathUtilities.tol8;
				boolean cx = (x[0] >= d4 - wp / 2) && (x[0] <= d4 + wp / 2);
				return cx && cy;
			}
		};

		double P = -100;

		/*
		 * distributed load on cap
		 */
		ArrayList<Node> left_force_nodes = MeshUtilities.seekNodes(mesh, steel_cap4);
		ArrayList<Node> right_force_nodes = MeshUtilities.seekNodes(mesh, steel_cap3);
		for (Node n : left_force_nodes) {
			n.setSortingCriteria(SortingCriteria.SORT_ON_X_COORDINATE);
		}
		for (Node n : right_force_nodes) {
			n.setSortingCriteria(SortingCriteria.SORT_ON_X_COORDINATE);
		}

		Collections.sort(left_force_nodes);
		Collections.sort(right_force_nodes);

		FunctionRnToR fy1 = new ConstantFunctionRnToR(P);
		FunctionRnToR fy2 = new ConstantFunctionRnToR(loadFactor * P);

		for (int i = 0; i < left_force_nodes.size() - 1; i++) {
			Node n1 = left_force_nodes.get(i);
			Node n2 = left_force_nodes.get(i + 1);
			Edge e = mesh.addEdge(new Edge(n1, n2));
			e.setNumberOfFields(2);
			Load l = new DistributedLoad(e, null, fy2);
			m.addLoad(l);
		}
		for (int i = 0; i < right_force_nodes.size() - 1; i++) {
			Node n1 = right_force_nodes.get(i);
			Node n2 = right_force_nodes.get(i + 1);
			Edge e = mesh.addEdge(new Edge(n1, n2));
			e.setNumberOfFields(2);
			Load l = new DistributedLoad(e, null, fy1);
			m.addLoad(l);
		}

		/*
		 * concentrated load on caps
		 */
		// Node n1 = MeshUtilities.seekNodeNearest(mesh, steel_cap4, d1, h +
		// hp);
		// Force P1 = new Force(n1, 0, -loadFactor * P);
		// m.addLoad(P1);
		//
		// Node n2 = MeshUtilities.seekNodeNearest(mesh, steel_cap3, d4, h +
		// hp);
		// Force P2 = new Force(n2, 0, -P);
		// m.addLoad(P2);

		return m;
	}

	@Override
	public MaterialModel createMaterial(Model m) {
		double E = 28.8e3;
		double nu = 0.18;
		double t = 156;
		State ss = State.PLANE_STRESS;
		MaterialModel mm = new StVenantKirchhoffMM(E, nu, t, 0, ss);
		return mm;
	}

	public MaterialModel createMaterial2(Model m) {
		// double E = 28.8e4;
		// double nu = 0.18;

		double E = 20e6;
		double nu = 0.3;

		double t = 156;
		State ss = State.PLANE_STRESS;
		MaterialModel mm = new StVenantKirchhoffMM(E, nu, t, 0, ss);
		return mm;
	}

	@Override
	public Model createModel() {

		Mesh mesh = createMesh();

		/*
		 * set the number of fields for structural problem
		 */
		mesh.setNumberOfFields(mesh.getDimension());

		Model m = new Model(mesh);

		MaterialModel mm = createMaterial(m);
		MaterialModel mm2 = createMaterial2(m);

		CohesiveLaw cl = createCohesiveLaw(m);

		// SpatialCriteria steel_cap_region = new SpatialCriteria() {
		//
		// double[] c1 = new double[] { 600.0, 0.0 };
		// double[] c2 = new double[] { 1119.0, 0.0 };
		// double[] c3 = new double[] { 722.0, 306.0 };
		// double[] c4 = new double[] { 203.0, 306.0 };
		// double r1 = 25.0;
		// double r2 = 25.0;
		//
		// @Override
		// public boolean isSatisfied(double... x) {
		// boolean is = false;
		//
		// /*
		// * check if in region of cap 1
		// */
		// is = is || (Math.sqrt(Math.pow(x[0] - c1[0], 2) + Math.pow(x[1] -
		// c1[1], 2)) < r2);
		// /*
		// * check if in region of cap 2
		// */
		// is = is || (Math.sqrt(Math.pow(x[0] - c2[0], 2) + Math.pow(x[1] -
		// c2[1], 2)) < r1);
		// /*
		// * check if in region of cap 3
		// */
		// is = is || (Math.sqrt(Math.pow(x[0] - c3[0], 2) + Math.pow(x[1] -
		// c3[1], 2)) < r1);
		// /*
		// * check if in region of cap 4
		// */
		// is = is || (Math.sqrt(Math.pow(x[0] - c4[0], 2) + Math.pow(x[1] -
		// c4[1], 2)) < r1);
		//
		// return is;
		// }
		// };

		SpatialCriteria steel_cap_region = new SpatialCriteria() {

			double[] c1 = new double[] { 600.0, 0.0 };
			double[] c2 = new double[] { 1119.0, 0.0 };
			double[] c3 = new double[] { 722.0, 306.0 };
			double[] c4 = new double[] { 203.0, 306.0 };
			double wp = 40;
			double hp = 20;

			@Override
			public boolean isSatisfied(double... x) {
				boolean is = false;

				/*
				 * check if in region of cap 1
				 */
				is = is
						|| ((x[0] <= c1[0] + wp / 2) && (x[0] >= c1[0] - wp / 2) && (x[1] <= c1[1]) && (x[1] >= c1[1]
								- hp));
				/*
				 * check if in region of cap 2
				 */
				is = is
						|| ((x[0] <= c2[0] + wp / 2) && (x[0] >= c2[0] - wp / 2) && (x[1] <= c2[1]) && (x[1] >= c2[1]
								- hp));
				/*
				 * check if in region of cap 3
				 */
				is = is
						|| ((x[0] <= c3[0] + wp / 2) && (x[0] >= c3[0] - wp / 2) && (x[1] >= c3[1]) && (x[1] <= c3[1]
								+ hp));
				/*
				 * check if in region of cap 4
				 */
				is = is
						|| ((x[0] <= c4[0] + wp / 2) && (x[0] >= c4[0] - wp / 2) && (x[1] >= c4[1]) && (x[1] <= c4[1]
								+ hp));

				return is;
			}
		};

		for (int i = 0; i < mesh.countFaces(); i++) {
			Face f = mesh.getFace(i);
			Element e;
			/*
			 * check if the face is in the region of steel cap
			 */
			double[] c = f.computeCenter();
			if (steel_cap_region.isSatisfied(c)) {
				/*
				 * enhance steel cap region with stronger material
				 */
				e = createElement(f, mm2, cl);
			} else {
				/*
				 * the other region is weaker
				 */
				e = createElement(f, mm, cl);
			}

			m.addElement(e);
		}

		m = createConditions(m);

		return m;
	}

	/**
	 * create the strategy
	 * 
	 * @param m
	 * @return
	 */
	protected abstract Strategy createStrategy(Model m);

	@Override
	public void addObservers(final Model m, final Strategy s, TimeStepObserver o) {
		addFundamentalObservers(m, s, o);
		addObserversForCrackProblem(m, s, o);

		s.addObserver(new StressObserver(m));

		XYSourceObserver xyo = new XYSourceObserver() {

			Node n1 = MeshUtilities.seekNodeNearest(m.getMesh(), 654.0, 0);
			Node n2 = MeshUtilities.seekNodeNearest(m.getMesh(), 668.0, 0);

			@Override
			public void fireTimeStepFinished() {
				double[] u1 = n1.getUHatLocal();
				double[] u2 = n2.getUHatLocal();
				double du = LinAlg.diffNorm2(u1.length, u1, u2);
				double lambda = s.getCurrentTimeStep();
				setData(du, lambda);
			}

		};

		s.addObserver(xyo);
		XYGraphicalObserver xygo = new XYGraphicalObserver(xyo, "lambda-distance plot");
		String fn = projectDir + '/' + projectName + '/' + meshFileName + ".ldcurve.txt";
		try {
			xygo.setOutputToFile(fn);
		} catch (IOException e) {
			e.printStackTrace();
			throw new Error("set output failed");
		}
		s.addObserver(xygo);
	}

	@Override
	public void run() {
		Model m = createModel();

		Analysis an = new Analysis(m);

		Strategy s = createStrategy(m);

		s.setEchoLevel(EchoLevelType.OUTPUT, EchoLevelType.AUX1);

		an.setStrategy(s);

		addObservers(m, s, (TimeStepObserver) s.addObserver(new TimeStepObserver()));

		an.run();

		// postProcess(m);

		Debugger.warn("Analysis completed");
	}

}
